output_version <- '2024-07-18'
source(paste0(wkdir,'/udf/regional_seas.R'))

## fishnet
fishnet_lake_chad_ntl_adm0_x_d01_shp <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01')
fishnet_adm0_only <- read_sf(dirname(fishnet_lake_chad_ntl_adm0_x_d01_shp),basename(fishnet_lake_chad_ntl_adm0_x_d01_shp))
st_crs(fishnet_adm0_only) <- "EPSG:4326"

##
## Groundwater depth
##

gwdepth_fn <- file.path(wkdir,'local/gis_data/Groundwater_Depth/DepthToGroundwater2021/Finalcombined_map/dtgw_20211.tif')
gwdepth_raw <- terra::rast(gwdepth_fn)
is.factor(gwdepth_raw)
factor_raster <- (as.factor(raster::raster(gwdepth_raw)))

#VS 0-7 - 254,972 obs
#S 7-25 - 167,426
#SM 25-50 - 164,138
#M 50-100 - 121,364
#D 100-250 - 144,664
#VD > 250 - 150,307 obs
labels0 <- c("VS 0-7", "S 7-25", "SM 25-50", "M 50-100","D 100-250","VD > 250")  # Define your labels

# Specify your categories; you can customize this based on your needs
rat <- levels(factor_raster)[[1]]
rat$dtgw_20211 <- labels0
levels(factor_raster) <- rat

require(exactextractr)

sum_cover <- function(x){
  list(x %>%
         group_by(value) %>%
         summarize(total_area = sum(coverage_area)) %>%
         mutate(sh_gw_depthcat = total_area/sum(total_area)))
  
}


##
## RADIUS BUFFER 0 (grid cell only) 50KM (FROM POINT), 100KM (FROM POINT)
##

library(doParallel)
cl <- makeCluster(12)
registerDoParallel(cl)

n_objectid_list=unique(fishnet_adm0_only$OBJECTID)
# n_objectid_list
results <- foreach(objectid_j=n_objectid_list,
                   .export=c(ls()),
                   .packages = c("dplyr", "sf","exactextractr","spatialEco")) %dopar% {
  # Modulus operation
  if(objectid_j %% 100==0) {
                       # Print on the screen some message
                       cat(paste0("iteration: ", objectid_j, "\n"))
  }
  
  fishnet_lake_chad_ntl_adm0_x_d01_shp_fn <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01.shp')
  fishnet_adm0_only_j = read_sf(fishnet_lake_chad_ntl_adm0_x_d01_shp_fn, query = sprintf('SELECT * from fishnet_lake_chad_ntl_adm0_x_d01 WHERE OBJECTID = %s',objectid_j))
  st_crs(fishnet_adm0_only_j) <- "EPSG:4326"
  
  if(dim(fishnet_adm0_only_j)[1]>0){
    ## 0km
    x <- exact_extract(factor_raster, fishnet_adm0_only_j, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover)
    
    #merge the list elements into a df
    gw_depth_share0 <- bind_rows(x, .id = "objectid") %>%
      mutate(buffer_km = 0)
    
    remove(x)
    # projected 50km
    fishnet_adm0_only_buf50km <- spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=50000)
    x <- exact_extract(factor_raster, fishnet_adm0_only_buf50km, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover)
    
    #merge the list elements into a df
    gw_depth_share50 <- bind_rows(x, .id = "objectid") %>%
      mutate(buffer_km = 50)
    
    # projected 50km
    fishnet_adm0_only_buf100km <- spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=100000)
    x <- exact_extract(factor_raster, fishnet_adm0_only_buf100km, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover)
    
    #merge the list elements into a df
    gw_depth_share100 <- bind_rows(x, .id = "objectid") %>%
      mutate(buffer_km = 100)
    
    out <- rbind(gw_depth_share0,gw_depth_share50,gw_depth_share100)
    
    # add objectid
    out$objectid <- objectid_j
    
    # clean-up
    rm(x,fishnet_adm0_only_buf100km,fishnet_adm0_only_buf50km)
    
    return(out)  
  }
    
}

parallel::stopCluster(cl)

system.time(gw_dt <- data.table::rbindlist(results))
gw_df <- gw_dt %>% as.data.frame()

# change name for 
names(rat)[2] <- "gw_depthcat"

gw_depth_out <- merge(gw_df,rat,by.x='value',by.y='ID',all.x=T) %>%
  dplyr::select(objectid,sh_gw_depthcat,gw_depthcat,buffer_km) %>%
  arrange(objectid,gw_depthcat,buffer_km)

gw_depth_out_fn <- file.path(wkdir,'proc_data',sprintf('GWDEPTH_share%s.dta',output_version))

# Adding the label value
attr(gw_depth_out$objectid, "label") <- "objectid of grid"
attr(gw_depth_out$sh_gw_depthcat, "label") <- 'Share of area covered by Groundwater Depth category : Bonsor and MacDonald 2011'
attr(gw_depth_out$gw_depthcat, "label") <- 'Groundwater Depth category : Bonsor and MacDonald 2011'
attr(gw_depth_out$buffer_km, "label") <- 'Buffer'
attr(gw_depth_out, "label") <- sprintf("sh_area_gw_depth_cat.R last run on %s",Sys.Date())
head(gw_depth_out)
summary(gw_depth_out)

## Save the data
haven::write_dta(gw_depth_out, gw_depth_out_fn)

# clean-up
rm(gw_dt)
rm(gw_df)
rm(results)